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Two-dimensional Boussinesq convection is studied numerically with very fine spatial resolutions 
up to 4096 2 . Our numerical study starts with a smooth asymmetric initial condition, which is 
chosen to make the flow field well confined in the computational domain until the blow-up time 
(T c ). Our study shows that the vorticity will blow up at a finite time with |w| mal ~ (T c — t)~ 1,61 
and \V6\ max ~ (T c -fr 3 - 58 . 

PACS numbers: 47.20.Cq, 47.27.Te, 47.27.Eq, 47.27.Jv 



Understanding whether smooth initial conditions in 
three-dimensional (3D) Euler equations can develop sin- 
gularities in a finite time is an important step in under- 
standing high- Reynolds- number hydrodynamics Q, 0- 
Two-dimensional (2D) Boussinesq convection is a sim- 
plified model of the 3D axisymmetric flows with swirl 
if main flow structures are well away from the symme- 
try axis ||. The computing requirements of the 2D 
Boussinesq simulations are significantly less than those 
of the 3D Euler equations. However, the spatial resolu- 
tions adopted have still been not fine enough to resolve 
the small scale structure. Therefore, the spectral meth- 
ods and the adaptive mesh methods H, IE Hi are 
mostly adopted. One central issue of this study is on 
whether there exist singularities in the 2D Boussinesq 
flows. If we denote T c the time of blowup, the min- 
imum criterion for the breakdown of smooth solutions 
in the 2D Boussinesq equations is: |w| maa ; ~ (T c — t) 
and \V9\ max ~ (T c - t)~ P , where a > 1 and (3 > 2 
0, • Adaptive mesh techniques 0. 0. 0] and the spec- 
tral method 3 are adopted to investigate this problem, 
which yield various conclusions, e.g., 0, observe no 
vorticity blow-up, while |3J,|5( only provide the marginal 
values (a — 1 and (3 = 2) although these studies predict 
vorticity singularities. 

Since the fast developed adaptive mesh methods are 
limited by the finite-order accuracy and mesh equality, a 
further "spectral" effort seems necessary to investigate 
this challenging problem. The spectral methods have 
been used intensively in 3D Euler simulations, but they 
are limited by the available computing capability: the 
finest resolutions used so far have been 2048 3 M The 
conclusions are affected by certain kind of symmetric 
assumptions introduced to increase the effective resolu- 
tions 0, ^3 • We can argue that the conclusion drawn by 



the 3D studies 0, E| is n °t necessarily more convincing 
than the axisymmetric assumption 4] where a resolu- 
tion of 1500 2 is adopted. In this study, spectral com- 
putations with extremely high resolution (from 2048 2 to 
4096 2 ) are carried out to investigate the blow-up issue 
for the 2D Boussinesq convection problem. Moreover, 
we try to maintain the flow structure well away from the 
axis at the time when solutions begin to blow up. This is 
done by choosing appropriate initial data. The govern- 
ing equations are solved by a fully de-alias pseudospectral 
Fourier methods with 8/9 phase-shifted scheme. The dig- 
ital filter is adopted to modify the Fourier coefficient to 
increase the stability of the numerical codes. The ma- 
chine accuracy of our computer with double precision is 
e = 10~ 16 « e~ 37 , and the modifying factor in the filter 
is <p(k) = e" 37 ^) 16 for k < N 12]. With these efforts, 
the vorticity blow-up is observed, which is in contrast 
to the conclusion drawn by the earlier spectral computa- 
tions 3 . 

The non-dimensionized 2D inviscid Boussinesq convec- 
tion equations can be written in the uj-tp formulation: 



e t + u • we = o, 

Wt + U' X7uj = - 
Aip = -oj, 



(1) 
(2) 
(3) 



where the gravitational constant is normalized to g = 
(0,-1), is the temperature, u = (u,v) the velocity, v 
the kinematic viscosity, uj — (0, 0, to) — V x u vorticity, 
and ip stream function. The simulation is carried out in 
the [0, 27r] 2 doubly-periodic domain. At first, we take the 
initial condition with unified zero vorticity and a cap-like 
contour of temperature with the following expression: 



0(x,y,O) = ( 



Ax — 3tt 



)6 1 (x,y)6 2 (x,y)[l-e 1 (x,y)], (4) 
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where if S(x, y) := 7r 2 — y 2 — (x — n) 2 is positive, 
then 6*i = exp (l — tt 2 /S(x, y)) , and zero otherwise; 
if s(y) := \y — 2ir\ /1.95-7T is less than 1, then 82 — 





exp (l — (1 — s(?/) 2 ) _1 ) , and zero otherwise. This ini- 
tial condition jtj is similar to except that a factor 
(Ax — 37r)/7r s used to break down the field symmetry 
with respect to x = tt. The main flow structure looks 
like a rising bubble which touches the y = and y = 2tt 
boundaries at t ~ 2.0. This will cause the singularity at 
the axis if we transform the 2D Boussinesq results back 
to the 3D axisymmetric flow [j] . To fix this problem, we 
compress the intermediate results at t = 1.2 obtained by 



using the initial condition to form a new initial data. 
More precisely, we let w(x, y, 0) = ui'(x, 2y— OAtt, 1.2) and 
d(x,y,0) = 6'(x,2y-QAw,1.2), for (x,y) 6 [0, 2tt] x [0, tt] 
(where 8' and u' are obtained by solving Q-Q) with a 
2048 2 grid), and zero otherwise. The new initial condi- 
tion only needs about 1/4 simulation time to reach T c 
comparing with the run in Q , and advantages of higher 
resolutions show up at early stages of simulations (Figs. 
Q). The time steps for the resolutions 1024 2 , 2048 2 and 



FIG. 3: The evolution of the T 2 , T4 and maximum 8 errors for three different resolutions (I024 2 : solid line, 2048 2 : 
dashline, 4096 2 : circle line). The errors are respectively defined as (T 2 (0) - T 2 (t))/T 2 (0), (T 4 (0) - T 4 (t))/T 4 (0) and 

\9max(0) — 8 m ax(t)\/\0m ax (0)\. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



4000 



&_ 1000 



500 




0.2 0.4 0.6 0.8 



FIG. 4: (a) and (b) are the time evolutions of maximum 
M and |V0| (1024 2 : solid line, 2048 2 : dashline, 4096 2 : circle 
line). It seems that 1024 2 and 2048 2 runs have almost iden- 
tical peak values from t = until t = 0.15, and the results 
of 2048 2 and 4096 2 are very close until t = 0.3. There are 
two reasons for these departures: 1) the filter removes more 
energy for lower resolution; 2) higher resolution can make a 
better fit for delta-like functions. 



4096 2 are 0.0004, 0.0002 and 0.0001, respectively, given 
by the CFL condition. 

The compressed bubble (Figs. IHa)(d)) will continue 
to rise, and the edge of the cap will roll up at later 
times (Figs. dc)(f)). At t ~ 0.9, the density and 
vorticity contours develop into the shape of "two asym- 
metric eyes." The contour revolutions on three differ- 
ent resolutions reveal similar phenomena at this point. 



On the other hand, combining the divergence-free condi- 
tion, the doubly-periodic condition and Eq. fll"|l. we can 
verify that T 2 (t) = 6 2 {x,y,i)dxdy and T 4 (<) = 

Jo ^ Jo* ^ 4 ( a; ' Vi t)dxdy are time independent. Figs. El (a) 
and (b) demonstrate that the global average quantities 
are well conserved within 1% error for all the three res- 
olutions used. The flow field develops into many small 
vortices after t = 0.8 and more and more energy is re- 
moved by the filter in the time evolution. Hence, it seems 
that simulations with higher resolution do not have more 
advantages here. 

However, the time evolution of the maximum values 
of 6, which is also time-independent, shows much better 
improvement for the 4096 2 run by comparing with lower 
resolution runs (Fig. 13(c)): the relative error for |0| maa; 
with the 4096 2 run is always lower than 10 -4 , while the 
corresponding error with the 1024 2 run is about 10~ 2 . 
Because the global maximum values of \u\ and |V0| are 
normally used as the indicators in the singularity analy- 
sis, the 4096 2 run may lead to a more accurate conclu- 
sion than that of the 1024 2 and 2048 2 runs. It should be 
pointed out that if resolutions higher than 4096 2 are em- 
ployed then the quadruple precision (e = 10 -32 ) instead 
of the double precision (e = 10~ 16 ) may have to be used; 
otherwise the results may be spoiled by the round-off er- 
rors (see [3 for more details). The maximum vorticity 
is initially located on the left edge of the rising bubble 
(Fig. 12a)), and later moved to the outer layer of the 
"eyes", (see, Figs, ^c) and|2Ia). This observation is in 
good agreement to the prediction in [4J. At t = 0.88, the 
\uj\ max locates at (5.31, 2.72) (see the "+" in Fig. 0(a)). 

Around t ~ 0.88, some very small vorticity structures 
begin to appear at the lower part of the smooth outer 
layer where the maximum \lu\ appears. The filter we 
adopted in the code will remove more and more energy 
from the system, so that the global average values like 
T 2 and T 4 are greatly affected (Figs. E](a)(b)). Conse- 
quently, after this critical time the numerical results in 
the 4096 2 run may not be accurate and reliable. Actu- 
ally, it is observed that the |w| ma:E and |V0| maa; expe- 
rience a drop-down after t ~ 0.9. This indicates that 
the filter prevents the blowup of the maximum vorticity 
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FIG. 5: The time evolutions of re-scaled |w|ma:c (square line), 
|V0|maa; (star line) and their least square fit to a straight line 
(dash line) for the (a) 2048 2 and (b) 4096 2 run. The blow-up 
times (To) are the points where x-axes are intersected by dash 
lines. The 1024 2 run has not been plotted because there is no 
blow-up signal until t = 1.0. 

from happening, which is similar to the viscosity effect 
in high Reynolds number simulations (e.g. 0, [l4^). 
Therefore, the sample maximum values after the drop- 
ping down should not be used in the singularity analysis. 
For the 1024 2 and 2048 2 runs, however, the drop-downs 
appear later than t ~ 0.88 and the blow-up time T c is 
also delayed. In Figs. we only provide time evolu- 
tions for |w| ma a; and |V#| ma:I ; in different resolutions be- 
fore they drop down. We re-plot these maximum value 
evolutions in Fig. [SJ It seems that the growth of |cj| ma:E 
and \V9\ max in the 4096 2 run (Fig. EJt>)) plausibly ter- 
minates in a finite time with |cj|„ la:E ~ (T c — t)~ ' , and 



crudely | V9\ max ~ (T c — t) 3 58 (in the later case we only 
adopt the values of \V6\ max after t — 0.75). It is no- 
ticed that our values (a,/3) = (1.61,3.58) are larger than 
the minimum blow-up values 0, Q and the existing re- 
sults We run our simulations until t = 2.0 and 
find that the main structure of the flow field stays away 
from y = and y = 2tt. The axisymmctric assumption 
adopted in the original 3D Euler equation model does not 
cause any axis singularity difficulty as discussed earlier 
with our new initial condition. 

We have not performed the simulation on grid finer 
than 4096 2 , but we can predict some results from the 
present computations. It is well known that when a 
delta function is approximated by a finite Fourier series, 
each doubling of the resolutions will cause doubling of 
the maximum value and 2™ +1 times the maximum val- 
ues of the n th order space derivative (see the analysis in 
Appendix A of ^2). ^ n the final state of our simula- 
tion (Fig. EI a)), the cut-line at y = 2.72 through the 
out layer of the "eye" looks very similar to a delta func- 
tion. Therefore, when finer and finer resolutions are used 
the peak vorticity and temperature gradient are getting 
larger and larger. Further 2D Boussinesq simulations 
with even higher resolution will support our singularity 
prediction although it seems impossible to accurately pre- 
dict the singularity time with the current computational 
scheme. In Fig. 0{b), the gradients of the curves become 
larger when the grids are refined. The 1024 2 result in- 
dicates no blow-up, and the singularity analysis for the 
2048 2 run (Fig. EJa)) shows that a = 1.28, = 2.32 and 
T c = 1.47. The results obtained by the three resolutions 
reveal that the values of T c may be smaller than 1.34, and 
the values of a and (3 may be larger than 1.61 and 3.58 if 
much fine resolutions (> 4096 2 ) are used. Although our 
simulations can not provide accurate value for T c , the 
existence of singularity in the 2D Boussinesq simulations 
is strongly supported by our numerical computations. 
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